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^ Abstract 

^ The decay of a passive scalar in a three-dimensional chaotic flow is studied using high-resolution 

numerical simulations. The (volume-preserving) flow considered is a three-dimensional extension of 
^1 the randomised alternating sine flow employed extensively in studies of mixing in two dimensions. 

lO It is used to show that theoretical predictions for two-dimensional flows with small diffusivity carry 

over to three dimensions even though the stretching properties differ significantly. The variance 

decay rate, scalar field structure, and time evolution of statistical moments confirm that there are 
two distinct regimes of scalar decay: a locally controlled regime, which applies when the domain 
size is comparable to the characteristic lengthscale of the velocity field, and a globally controlled 
regime, which when applies when the domain is larger. Asymptotic predictions for the variance 

(N 

^ decay rate in both regimes show excellent agreement with the numerical results. Consideration of 



U 



as 
in 



both the forward flow and its time reverse makes it possible to compare the scalar evolution in 

flows with one or two expanding directions; simulations confirm the theoretical prediction that the 
decay rate of the scalar is the same in both flows, despite the very different scalar field structures. 
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I. INTRODUCTION 



In studies of fluid mixing, much attention has been devoted to chaotic-advection flows, 
characterised by relatively simple (smooth and divergence-free) velocity flelds but complex, 
chaotic particle trajectories. Flows of this type appear in several contexts: Stokes flows [1], 
turbulent flows with large Schmidt number (in the so-called Batchclor regime [2]), geophys- 
ical flows [3, 4], and elastic turbulence [5]. One aspect is of particular interest: the decay, 
through the combined effect of advection and diffusion, of the concentration of passive scalars 
released in such flows. In bounded domains this decay is exponentially fast; one of the main 
problems is then to relate the corresponding decay rate to the flow characteristics and to 
the diffusivity. 

For small diffusivity k — >■ 0, a clear distinction emerges between flows that are everywhere 
chaotic (e.g. uniformly hyperbolic), and flows with regular regions (e.g. containing elliptic 
islands [6] or incorporating no-slip boundary conditions [7, 8]): the decay rate of passive 
scalars tends to a non-zero limit as k — > in the flrst case, whereas it tends to zero in 
the second case. In this paper, we concentrate on the first case and, more specifically, on 
stochastic models of chaotic flows. These models, in which the complex time dependence of 
realistic flows is represented by random processes [9, 10] , have been the subject of sustained 
research from the mid 1990s onward [11]. Remarkably, several theoretical results are now 
available which relate properties of the scalar decay, notably the decay rate, to flow properties 
such as stretching statistics [12-19]. These results, which we review in Sec. II below, have 
been confirmed by numerical simulations [18, 19]. Confirmation, however, has been limited 
to two-dimensional fiows. The main aim of the present paper is to test the theoretical results 
against numerical simulations of three-dimensional flows. 

The passage from two to three dimensions adds signiflcantly complexity to the problem. 
Numerically, the high resolutions needed to capture the flne scalar scales generated for small 
diffusivitics make the 3D case much more computationally demanding than 2D. Physically, 
the dynamics of scalars in 3D flows are also much richer than in 2D flows: incomprcssibility 
provides a strong constraint in 2D, leading in particular to equal and opposite (stretching) 
Lyapunov exponents; in 3D incomprcssibility represents a looser constraint, imposing only 
that the sum of the three Lyapunov exponents vanish. The existence of three distinct 
Lyapunov exponents also has a profound influence on the concentration fleld. Depending 
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on whether the intermediate Lyapunov exponent is positive or negative, the concentration 
field is dominated either by quasi-two-dimensional structures ("pancakes") or by quasi- 
one-dimensional structures ("needles"). Interestingly, the most recent theoretical results 
about the scalar decay rate [18, 19] are insensitive to this difference in structure and make 
predictions that are independent of the sign of the intermediate Lyapunov exponent. Our 
simulations confirm the validity of this conclusion. 

A well-known difference between mixing in 2D and 3D is that chaotic particle trajecto- 
ries are only possible in 2D fiows if they are time dependent, whereas they exist in time- 
independent 3D flows. Time-independent 3D flows, however, are similar to one-degree-of- 
freedom Hamiltonian maps [20], and do not lead to global chaos but rather to a mixture of 
chaotic and (possibly small) regular regions, with the latter controlling scalar decay in the 
long-time limit. The flows with random time dependence that we consider do not have this 
type of intricate structure: they are statistically homogeneous and globally mixing. This 
is one of the aspects which differentiate this paper form earlier work on mixing in 3D, in 
particular the numerical simulations of time-independent ffows in [21] and [22]. 

Qualitatively, the mixing properties of sufficiently chaotic ffows are independent of the 
flow details. It is therefore advantageous to devise a random-flow model that is easy to 
implement and analyse. In 2D, this role is played by the randomised sine map [23, 24] 
which has been used extensively in studies of passive and reactive scalars. In this paper, we 
propose a straightforward 3D extension of the sine map. With the parameters chosen, this 
map has two positive and one negative Lyapunov exponent, which yield concentration flelds 
with pancake-like structure. We use the associated inverse map as a model for ffows with 
two negative and one positive Lyapunov exponent which generate needle-like concentration 
flelds. The extension to 3D is not uniquely deflned. Other versions of the 3D map are 
less mixing in that they do not stretch certain flxed directions. We briefly comment on the 
impact this property has on the scalar decay. 

The plan of the paper is as follows. In Sec. II, we review the theoretical predictions 
for the decay rate of the scalar concentration. Following [19], we emphasise the existence 
of two regimes of scalar decay: a locally controlled regime relevant to flows whose typical 
scale is comparable to the size of the (periodic) domain considered, and a globally controlled 
regime relevant to flows of smaller scale (see also [18]). The three-dimensional sine map is 
introduced in Sec. III. There we examine the stretching properties of the map and describe 
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the numerical procedure employed for simulations of the scalar evolution. Sees. IV and V 
present results of simulations carried out in different domain sizes. This makes it possible 
to explore both locally and globally controlled regimes and contrast differences in the decay 
rate as well as in the structure of the concentration field. We also compare concentration 
fields obtained with the map and its inverse. The paper concludes with a discussion in 
Sec. VI. 

II. DECAY RATE 

The evolution of the concentration C{x, t) of a passive scalar released in a flow v{x, t) is 
governed by the advection-diffusion equation 

dtC + v-VC^ kV^C, (1) 

with initial conditions C{x,0) — Co{x). When non-dimensionalised using characteristic 
length and velocity scales, this problem retains its form, with k now representing the inverse 
Peclet number rather than the diffusivity. This non-dimensional interpretation will be used 
in what follows. 

We focus on incompressible fiows, V • = 0, generated by specific random processes. 
These flows are assumed to have homogeneous and stationary statistics and to be spatially 
smooth, that is, \\v{x) — v{x')\\ oc ||a; — x'\\ as \\x — x'\\ — >■ 0. Periodic boundary conditions 
are applied to (1). Since (1) is left unchanged when a constant is added to C, there is no 
loss of generality in assuming that the concentration fleld averages to zero: 

J C{x,t)dx^ J Co{yi)dx^O. (2) 

By analogy with the flnite-dimensional case the concentration C{x,t) can be expected to 
decay exponentially in the long-time limit: 

C{x,t) ^ D{x,t)e~^^ as ^►oo (3) 

for almost all realisations of v{x,t), Here the deterministic decay rate A is best interpreted 
as (minus) the Lyapunov exponent of the linear system (1), and D{x,t) is a stationary 
function termed strange eigenmode [23] (see [25] for rigorous results). We emphasise that A 
is the Lyapunov exponent of the inflnite-dimensional random system (1) and should not be 
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confused with the Lyapunov exponent of the three-dimensional hnear system governing the 
separation of nearby particles in the velocity field v. We also note that is not the only 
useful time scale characterising the scalar decay: a dissipation time scale can for instance 
be defined as the time for some norm of C to be reduced from its initial value by a given 
fraction [26]; this, however, characterises the early-time behaviour rather than the long-time 
behaviour on which we focus. 

A striking feature of the decay rate is that it tends to a non-zero value in the limit of 
zero diffusivity (infinite Peclet number), limK_^o A = Aq 7^ 0, provided the fiow is sufficiently 
mixing. This raises the question of the relationship between the asymptotic decay rate, Aq, 
and the statistical properties of the flow, v(x,t). 

In a series of papers [11-14], Aq has been related to the large-deviation statistics of the 
finite-time Lyapunov exponents of v{x,t). These are defined as 

h = t-Hog\\y\\, (4) 

where y denotes the separation between nearby trajectories, and are distributed according 
to a pdf, t) say. Aq may be expressed in terms of the rate function 

g{h) ^ - Urn log p{h;t) (5) 

t— >oo 

associated with this pdf. Alternatively, it may be expressed in terms of the generalised 
Lyapunov exponent 

£(q) = lim rMogE||2/||«, (6) 

where E denotes ensemble average. These formulations are equivalent since g{q) and i{q) 
are related by a Legendre transform [27, 28]: 

i{q)= sup {qh-g{h)). (7) 

h 

Results of this type are not general, however. In particular, they do not apply to fiows 
whose typical scale is much smaller than the domain size [15, 16, 29, 30]. In fact there are 
two different regimes of scalar decay [17-19]. Flows whose spatial scale is comparable to the 
domain size are in a locally controlled regime, in which Aq depends only on the stretching 
statistics of v{x,t) and is given in terms of g{h) or i{q). Flows with smaller scale, on the 
other hand, are in a globally controlled regime, in which Aq depends on global fiow features 
including the domain size. 
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In [19], the two regimes are related to two different parts of the spectrum of the (de- 
terministic) hnear operator governing the evolution of the covariance EC{x,t)C{x' ,t) in 
Kraichan-Kazantsev and renewing flows [31]. The smallest non-zero eigenvalue determines 

7, the decay rate of the ensemble- averaged variance: 

7 = - lim log E / C^{x, t) dx. (8) 

t— >oo J 

In the locally controlled regime, 7 belongs to the continuous spectrum of the operator 
obtained for k — 0, while in the globally controlled regime, 7 represents a discrete eigenvalue. 

Because of this difference in the eigenvalues, the decay rates are qualitatively different in 
the two regimes. In the locally controlled regime the asymptotic form of 7 for «; ^ is 

7iocai = -£{q.) + + 0(1/ log' k), (9) 

log K 

where is the minimum of i{q) [so that i{q^) — —g{0), see (7)]. In the globally controlled 
regime, no such simple expression exists, and the leading-order term can be obtained only by 
solving the eigenvalue problem for the covariance. However, in the limit of flow scale much 
smaller than domain size, a useful estimate is provided by homogenisation theory, which 
approximates (1) by a diffusion equation with an effective diffusivity Kes [32], whence 

STrVcff 

7global ~ 7hom = — p , (lUj 

where L is the largest of the domain periods. 

Note that (9) and (10) describe the decay rate of the ensemble- averaged variance. The 
decay rate of a single realisation is given by 

7 = - lim log / C'{x, t) dx = 2A (11) 

on using (3). In general, 7 7^ 7, the difference marking the intermittency of the scalar 
decay [33]. For the flows considered here, differences between 7 and 7 are neghgible (see 
Sec. IV). Indeed we will confirm that (9) can be used to predict the variance decay in single 
realisations. (The equality of 7 and 7 in the homogenisation limit (10) is clear since the 
spatially homogenised equation is deterministic.) 

The theory summarised above holds in any dimension [19]. Until now, numerical verifi- 
cation (of (9) and (10) in particular) has been restricted to two dimensions. In this paper 
we assess the applicability to three-dimensional fiows. 
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The additional dimension introduces interesting complexities. In two dimensions, incom- 
pressibility strongly constrains the distribution of finite-time Lyapunov exponents, leading 
to the symmetry property i{q) = i{—q — 2). As a result, = —1 and the leading-order of 
(9), which in general reads 7 ~ —('{q*) = 5'(0), can be alternatively computed as 7 ~ —(,{—1) 
(it is in this form that the decay rate appears in some of the theories [12, 13, 18]). In di- 
mensions 0? > 2, by contrast, there is no such symmetry property; rather what hold are the 
relationships [28] 

=£(-g-(i), and hence g~ (h) = g{-h) - dh (12) 

between the stretching statistics of a fiow and those of the time-reversed fiow (denoted 
by ^). A consequence is that (9) predicts that the decay rate of a scalar is the same in 
a flow and in the time-reversed flow. This is a remarkable prediction for c? = 3 because 
many features of the stretching in a flow and its time-reverse are different. In particular, the 
intermediate Lyapunov exponent changes sign; thus, if a flow has one contracting and two 
expanding directions, the time-reverse has one expanding and two contracting directions. 
As a result, the scalar flelds have very different structures in the two flows, with typical 
concentration isosurfaces taking the shape of (nearly two-dimensional) pancakes in one case 
and the shape of (nearly one-dimensional) needles in the other. Nonetheless, the scalar 
concentration decays at the same rate, as we verify below [34]. 

III. THREE-DIMENSIONAL SINE MAP 
A. Formulation 

The randomised two-dimensional sine map [23, 24] has become a standard tool for the 
study of mixing by spatially smooth flows. It is given by 

Xn+i ^Xn + a sin(y„ + 0„) , ^yn + b sin(x„+i + ipn), (13) 

where a and b are constant parameters and 0„ and ipn are independent random angles 
uniformly distributed in [0, 27r]. The domain considered is periodic so the right-hand sides 
in (13) are taken modulo 27r. This map provides the positions at times nr, n — 1,2 ■ ■ ■ oi 
particles advected by a succession of shear flows alternating in direction. Its application to 
(1) is straightforward: advection of a scalar held by (13) can be implemented very efficiently 
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using a lattice representation; diffusion for small k can be incorporated as smoothing by the 
heat kernel oc exp(— (x^ + y'^) / (Akt)) (or rather the periodised version thereof). 

In what follows, we use a straightforward three-dimensional generahsation of (13), namely 



Xn+l =Xn + a Sm{yn + 4>n), Vn+l = Un + h Sni{Zn + ^Jn) , = Zn + C sin(a;„+i + (fn) , (14) 

where a, b and c are constant, and 0„, ipn and are independent uniformly distributed 
random phases in [0,27r]. This map, which corresponds to successive applications of shear 
flows in the x-, y- and 2;-directions, preserves volume and has homogeneous statistics. As in 
two dimensions, these shear flows are simple enough that they can be integrated explicitly 
to yield (14). The stretching statistics of the map are controlled by the Jacobian matrix, 
given at the origin by 



^ 1 acos0„ 

lb cos ■0n 

^ c cos(a sin (p^ + <Pn) clc cos 0„ cos(a sin 0„ + (^„) 1 ^ 



(15) 



The three Lyapunov exponents are computed numerically from the singular values of 
AnAn-i ■ ■ ■ Aq. In what follows, we fix a = 6 = c = tt. In this case, the Lyapunov ex- 
ponents arc found to be 1.04, 0.58 and —1.62. Thus the map (14) has one contracting and 
two expanding directions. 
The inverse map, 

Zn+l^ Zn+CSm{Xn + iPn): Vn+l ^ Vn+b SiuiZn+l+l/jn) : ^^n+l = 2;„-Fa Sin(y„+i-F0„) , (16) 

corresponds to the time-reversed flow. It has Jacobian matrix A~^, and ior a — b — c — tt 
Lyapunov exponents 1.62, —0.58 and —1.04, hence one expanding and two contracting 
directions. 

Three-dimensional generalisations of (13) arc not uniquely determined. Maps diflcring 
from (14) and (16) via permutations of the variables can be constructed. However, they 
may not be transitive, in the sense that their Jacobian matrices have invariant directions 
independent of the random phases. Such maps may not be sufficiently mixing for predictions 
such as (9) to hold, even if the largest Lyapunov exponent is positive. We return to this 
point in Sec. VI. 
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FIG. 1. Generalised Lyapunov exponents for the map (14) with a = b = c = n. The results of two 
Monte Carlo computations (circles and squares) are shown together with the best fit by a parabola. 



B. Theoretical predictions 

We now consider the predictions (9)-(10) for the three-dimensional map (14) and its 
inverse (16). We choose the parameters a = b = c = tt and take r = 1. 

In order to find the coefficients in the asymptotic formula (9), we use a Monte Carlo 
approach to compute i{q) numerically in the vicinity of its minimum, q^. Because the ex- 
pectation, E 111/ 11^, is dominated by rare realisations, the sampling required is rather delicate; 
we therefore use importance sampling, namely the random resampling method described in 
[28]. In particular, we use an ensemble of 5 x 10^ reahsations to estimate i{q) over a range 
bracketing q^.., q G [—2,-1.4]; £(g*) and (!-"{q^) follow from a parabolic fit. Some numerical 
results are shown in Figure 1. Rough error bars are obtained by repeating this procedure 
10 times. This provides the estimates £(g*) = 1.022 ± 0.002 and f{q^) = 0.93 ± 0.09. 
Introducing these values into (9) then yields the estimate 



, 18.36(±1.78) , , 

7,oeai = 1.022(±0.002) + -\ >-. (17) 

log K 



for the variance decay rate in the locally controlled regime. This is tested in Sec. IV. 

The corresponding prediction for the inverse map (16) is exactly the same as that for the 
direct map (14). Indeed, Eq. (12) implies that (■~{q) is the reflection of l{q) about the line 
q — —3/2. Therefore, (-"{q*) — i{q*), {^~)"{q*) — ^"{q*) and hence 7iocai is unchanged. 

To examine both locally and globally controlled regimes of scalar decay, we consider 



the scalar evolution for the forward map in a periodic cube of (total) length 2t:P with 
P — 1,2- ■ ■] that is, the flow has unit cells of length 27r but the right-hand sides of (14) are 
taken modulo 27rP. The asymptotic result (9) is expected to hold for small values of P, while 
the homogenisation approximation (10) should be valid for P ^ 1. The effective diffusivity 
Keflf appearing in (10), is easily found by recalling that the variance of one of the coordinates, 
Xn say, satisfies Ex^ — 2Kesn; computing the variance from (14) gives Kes — Hence 

7global ~ 7hom = (18) 

for large P. 

We wish to estimate the value of P at which the transition from local to global control 
occurs. The predicted decay rate 7 is determined by whichever of the local and global decay 
rates is smaller. Comparing (17) with (18) we therefore expect local control for P = 1,2 
and global control for P > 3. 



C. Numerical procedure 

Numerical simulations are performed on a regular grid of A^^ points. The advection uses 
the 3-D sine map (14), while for numerical efficiency the diffusion is applied in spectral space 
using a Fast Fourier Transform. We have confirmed that an explicit finite-difference step 
yields essentially indistinguishable results. The code has been tested by reproducing results 
of [19] for the 2-D sine map (13). 

Two sets of initial conditions are considered. In most simulations we use 

C{x, 0) = sm{x/P) sin(y/P) sm{z/P). (19) 
In Sec. V, we also use the initial condition 

C{x, 0) = sm{z/P) (20) 

corresponding to one of the gravest modes of the diffusion operator in the domain. 

The simulations need to resolve spatial scales in the range between the box size, Lq, 
and the tracer microscale, L^,. The latter scale is estimated by matching the diffusive and 
stretching timescales, yielding = \/kJ\, where A is the largest Lyapunov exponent of v. 
Thus we require > Ax, where Ax = Lq/N is the grid spacing, or 

«>(|)^A. (21) 
10 



In the simulations of Sees. IV and Y, N — 1024. This imphes a critical value of the diffusivity, 



For P > 1, the numerical simulations are performed on the A^^ grid but with replicas 
of the velocity field. Taking P > 1 amounts to a redefinition of Lq. For simplicity, however, 
it is convenient to keep Lq = 27r fixed and rescale the equations. This allows us to reuse the 
code described above with minimal changes. More precisely, we use the modified map 



with the scahng k, i-> P'^k, and the right-hand sides taken modulo 27r. Note that the Jacobian 
is unchanged. 

IV. LOCAL CONTROL: P = 1 

In this section, we focus on the locally controlled regime and take P = 1. Wc report on 
simulations of the forward map (14) for several values of k in the range 10~^-10~^, and of 
the inverse map (16) for k — 10~^. 

The evolution of the variance 



in these simulations is shown in Fig. 2(a). As with the two-dimensional sine map (13), the 
variance decays exponentially after a brief initial transient, in agreement with the strange- 
eigenmode prediction (3). These results are insensitive to resolution: simulations at the 
lower resolution N — 512 (not shown) are almost identical. Also, there is little variability 
between different realisations. For an ensemble of 8 realisations evaluated at N — 512 (not 
shown), the decay rate of the variance for k, — 10~^ varies by about 2%. This difference 
is of the same order of magnitude as the difference between the decay rates for N = 512 
and = 1024. Wc therefore conclude that the intermittcncy is negligible and henceforth 
restrict attention to the single-realisation decay rate 7 which we identify with 7. 

The forward and inverse maps lead to almost identical decay rates. Specifically, for 
K — 10"^, the decay rates for the forward and inverse maps are 1.22 and 1.21, respectively. 
This supports the claim made in Sec. Ill B that their decay rates are equal despite differences 
in their stretching properties. 



«e = 0(10-5). 



Xn+1 ^Xn+ — sin(P?/„ + (t)n), Vn+l ^ Vn + p sin{PZn + ^n), 

Zn+1 ^ Zn+ — sin(Pa;„+i + cpn), 



(22) 




(23) 
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(a) 



(b) 




10 20 30 10"^ 10"= 10-* 10"^ 10"^ 10"' 



n K 

FIG. 2. Variance decay for P = 1. (a): variance versus n for the forward map with k = 
10~^, 10~^, 10~^, 10~^, 10~^ (sohd Unes, from top to bottom), and for the inverse map with k = 
10~^ (dashed Une). (b): decay rate 7 versus k for the forward map. The decay rates, obtained 

from a least-squares fit over n € [5, 30] (symbols, with error bars corresponding to the standard 
error), are compared with the prediction 7iocai in (17) (solid line, with dashed lines indicating the 
error estimates.) 



The dependence of 7 on k for the forward map is shown in Fig. 2(b). For each value of 
K, 7 is obtained from a least-squares fit over the time interval, n e [5, 30]. There is excellent 
agreement with the asymptotic estimate (17) for 7iocai- This confirms the applicability of 
the asymptotic result (9) to the locally controlled regime in more than two dimensions. The 
agreement deteriorates rather rapidly for larger k, but this is to be expected since the theory 
assumes k — >■ and neglects terms that are o(l/ log^ k). 

The spatial structure of the decaying scalar is illustrated by Figs. 3 and 4, which show 
volume-rendered concentration fields for the forward and inverse maps, respectively. The 
plots use the normalised concentration variable 

= (24) 

(7 

where a is the the standard deviation defined by (23). This provides a numerical approxi- 
mation to the strange eigenmode structure D{x,t) defined in (3). To facilitate comparison, 
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(a) (b) 




FIG. 3. (Colour online) Concentration field generated by the forward map (14) with K = 10-^: (a) 
n = 2; (b) n = 4; (c) n = 6. Note the presence of filamentary sheets caused by expansion in two 
directions and contraction in one direction. 



(a) (b) (c) 




FIG. 4. (Colour online) Concentration field generated by the inverse map (16) with k = 10~ : 
(a) n = 2; (b) n = 4; (c) n = 6. Note the presence of thin "needles" caused by expansion in one 
direction and contraction in two directions. 

the colour scale is renormalised in each panel, i.e., it spans the respective minimum and 
maximum values. 

For the forward map, there is a strong resemblance to the two-dimensional case: large- 
scale filaments dominate, though in the present case they are manifested as two-dimensional 
surfaces. This, of course, results from the presence of two expanding and one contracting 
directions. 
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For the inverse map, the picture is quite different. In place of the large-scale filamentary 
sheets, there are now thin "needles", which arise from the presence of two contracting 
directions. This is most obvious for n = 4. Note that these needles are approximately parallel 
to the (x, y)-plane: the final advection step of each iteration consists of a purely horizontal 
shear (sec (16)) which creates elongated structures in the {x, 7/)-plane. The contrast between 
Figs. 3 and 4 underscores the profound effect that differences in the stretching properties of 
the forward and inverse maps — in particular, the number of positive Lyapunov exponents 
— have on the structure of the concentration field. It is therefore remarkable that the 
variance decay rate is insensitive to this. The physical explanation for this result is that the 
variance for large n is controlled by isolated pockets of extreme concentration values; these 
pockets, which can be readily identified in Figs. 3c and 4c and correspond to rare events of 
low stretching [19], have the same statistics in the forward and inverse maps (cf. Sec. II). 

V. TRANSITION TO GLOBAL CONTROL: P > 1 

In this section we examine the transition between the locally controlled and globally 
controlled regimes by analysing the scalar decay in domains of size 2nP, P > 1. We focus 
on the forward map (14) and k — 10~^ and = 1024. Simulations have been carried 
out with the two initial conditions (19) and (20). 

A. Variance decay and concentration fields 

Fig. 5 summarises the variance decay results by showing the decay rate 7 as a function 
of P. The numerical results are compared with the asymptotic predictions (9) and (10) for 
the locally and globally controlled regimes. The existence of the two distinct regimes is clear 
from the figure, with the transition taking place around P = 3. The structure of the scalar 
field for P = 2 (see Fig. 6(a)), which is very similar to the structure for P = 1, suggests that 
P = 2 is in the locally controlled regime. This is consistent with the argument outlined in 
Sec. IIIB which predicts that P = 3 corresponds to the smallest domain for global control. 

The decay rates found for P > 3 with the initial conditions (19) are much larger than 
the approximation to 7giobai provided by homogenisation theory in (10). We attribute this 
discrepancy to finite-time or symmetry effects: in the homogenisation scenario, the decaying 

14 




FIG. 5. Decay rate 7 versus P for k = 10~^. The numerical results obtained with the initial condi- 
tions (19) [crosses, +] and (20) [asterisks, *] are compared with the prediction from homogenisation 
theory 7giobai (long dashes) and with the prediction Tiocai for the locally controlled regime (solid 
line) . 

concentration field has a structure close to that of one of the gravest modes of the diffusion 
operator, that is sin(a;/P), sin(y/P) or sin(2;/P), depending on details of the initial condi- 
tions. However, in our simulations with the initial conditions (19), the concentration fields 
are very different (see Figs. 6(b)-(d)). This can be understood as a finite-time effect. As 
the diffusive approximation of homogenisation theory indicates, the relative amplitudes of 
the Fourier modes depend on the initial conditions; therefore the high modes can dominate 
for some finite time if their initial amplitudes are substantially larger than the amplitudes 
of the gravest modes. This appears to be the case with the initial conditions (19), which 
would explain why the prediction (10), although valid for n ^ 1, is not realised in our 
simulations (which are limited to n < 30). Furthermore, for P > 2 even, the dynamics 
preserves the symmetries C{x -\- Ptt, z) = —C{x, y, z), C{x,y + Ptt, z) = —C{x, y, z) and 
C{x,y, z + Pn) = —C{x,y,z) of the initial conditions (19) (see Figs. 6(a),(c),(d)); however, 
these symmetries are not compatible with the gravest-mode structure. 

To confirm that the discrepancy between simulated and predicted decay rates results from 
finite-time effects and the special initial conditions, we have repeated the simulations using 
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(a) (b) 




FIG. 6. (Colour online) Concentration fields for P > 2, n = 6 and k = W ^ with initial conditions 
(19): (a) P = 2; (b) P = 3; (c) P = 4; (d) P = 8. Note the transition from a filamentary to a 
periodic appearance as P is increased. 



the gravest mode (20) as initial condition. In this case, the observed decay rates closely 
match those predicted by homogenisation theory. The scalar fields, displayed in Fig. 7, have 
the expected structure of a distorted gravest mode for P > 3, with the distortions reducing 
as P increases. Comparison between the scalar fields for P = 2 and P = 3 further supports 
our claim of local control for P = 2. 
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(a) 



(b) 



(c) 




FIG. 7. (Colour online) Concentration fields for k = 10~^ with initial conditions (20): (a) P = 2; 
(b) P = 3; (b) P = 4. As P increases, the structure of the concentration approaches the structure 
of one of the gravest modes of the diffusion operator, here s'mz. 

B. Statistical moments and pdfs 

An alternative view of the transition from local to global control is provided by statistical 
moments and pdfs of the concentration field. Normalised moments are defined by 



According to the strange-eigenmode expression (3), they should be stationary random func- 
tions. This was confirmed in numerical simulations of the 2-D sine map, which indicate that 
the mq are approximately constant [16]. This is at odds with the predictions of [13] that the 
rriq depend on q and increase exponentially with time. 

Fig. 8 plots TUq against n for various values of P and g = 4, 6, 8, 10. Focusing on P = 1, 
the rriq increase rapidly (approximately exponentially) on short timescales before levelling 
off, in close analogy to the two-dimensional results (see Fig. 3 of [16]). After the initial 
transient period, the rriq appear roughly independent of n, consistent with (3) and the 
existence of a statistical equilibrium established by stretching and diffusion. For P = 2 the 
picture is similar, though the initial growth of the ruq is smaller. For P = 4 and P = 8, 
by contrast, there is hardly any growth at all. This is consistent with the tracer behaviour 
being controlled by an effective diffusion for large P. 

With regards to the concentration pdfs, Fereday and Haynes [16] argued that the tails 




(25) 



17 
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10 20 30 10 20 30 

n n 

FIG. 8. Time evolution of the moments, niq vs n, for the forward map with K = 10-^: (a) P = 1; 
(b) P = 2; (c) P = 4; (d) P = 8. 

should scale like lO]"^ in the limit k — >■ 0, where 

/3~l + 2^(0)/7. (26) 

This implies that ^ ~ 3 in the locally controlled regime, since 7 ~ 7iocai ~ ~-^(?*) — 5'(0); on 
the other hand, ^ > 3 in the globally controlled regime, since 7giobai < Tiocai- This prediction 
was verified in the 2D case [16, 19]. Here we examine this prediction for the forward 3D 
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FIG. 9. p{9) vs 9 for the runs analysed in Fig. 8: (a) P = 1; (b) P = 2; (c) P = 4; (d) P = 8. 
The sohd Unes correspond to iterations n = 5, 10,20,30; the reference slope is plotted with 
a dotted line. In the global regime (panels c and d), the tails are much shorter. The pdfs are 
obtained by binning 9 G [—10, 10] with 200 bins of uniform width. 

map (14). 

Fig. 9 shows log-log plots of the pdf p{9) for the simulations analysed in Fig. 8. The 
solid lines correspond to pdfs evaluated at n = 5, 10,20,30, the dotted lines to |^|~^. We 
expect the spectral slope of the tails to change as P increases and we move from a local to 
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FIG. 10. (Colour online) Concentration field generated by the map (31). n = 3 and k = 10 ^. 

global regime. For P = 1 the slope is slightly steeper than the theoretical prediction, but 
the algebraic decay is clear enough. The discrepancy may be related to the finiteness of k. 
There is similar behaviour for P = 2. For P = 4 and P = 8, however, the prediction of 
algebraic tails breaks down. Here the pdfs decays rapidly for |^| > 1. This is consistent 
with the concentration field being controlled by an effective diffusion in the homogenisation 
regime P ^ 1. 



VI. DISCUSSION 

The simulations reported in this paper detail the structure of passive scalars decaying un- 
der the combined action of chaotic advection and (molecular) diffusion in 3D. They confirm 
that the strange-eigenmode behaviour (3), which has been well documented in 2D, carries 
over to 3D. By contrast with the 2D case, the strange eigenmodes are found to take three 
different forms, not only depending on whether the decay is locally or globally controlled, 
but also, in the former case, on the number of positive Lyapunov exponents. Support for 
this comes from statistical moments, concentration pdfs and fields, and the variance decay. 

The simulations confirm that the theoretical predictions of [19] for the decay rate of the 
concentration variance are valid in 3D as well as in 2D. In particular, the variance decay 
rate for — )• in the locally controlled regime is always given by g{0) = —i{q*) (Fig. 2), 
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the decay rate of the probabihty that a hne element experience no stretching in the time 
interval [0,t]. As a result, the variance decays at the same rate for the sine map (14) and 
its inverse (14), despite the very different structures of the scalar fields. The prediction for 
the variance decay rate in the globally controlled regime (9), which is obtained by applying 
homogenisation theory and assuming a large ratio of box size to characteristic scales, has 
also been well verified (Fig. 5). The transition from local to global control occurs when 
this ratio takes the value P = 3, in agreement with an estimate obtained by equating the 
predictions for the decay rates. 

The identical decay rates for the forward and inverse map may seem surprising. However, 
this can be explained straightforwardly. Physically, the variance decay is controlled by a few 
small fluid blobs that remain unstretched for long times. The fraction of fluid occupied by 
these anomalous regions decreases exponentially at a rate that is identical for the forward and 
inverse maps and unrelated to the number of positive Lyapunov exponents. Mathematically, 
the equality of the decay rates can be established directly, without resorting to the explicit 
expression (9) for the decay rate. To see this, let us denote by Tn the action of a volume- 
preserving map such as (14) on the concentration field, Cn = C{-,nT). Immediately after 
the nth advection step 

Cn+l{Xn+l) = {TnCn) (a^n+l) = Cn{Xn) (27) 

while after the nth diffusion step (i.e. completion of the full nth step) 

Cn+l{Xn+l) = {T^TnCn) {Xn+l), (28) 

where T) represents the effect of diffusive smoothing, Using volume preservation, it can be 
checked that the adjoint (in L2) of Tn is — (i.e. Tn is unitary), while — V. Now, 
after n + 1 steps of the forward map, the concentration field is 

Cn+l = VTnVTn-1 ' ' ' VToCq, (29) 

while after n + 1 steps of the inverse map it is 

Cn+i^V7^VTl,---VT^Co. (30) 

Comparing the operators on the right-hand side of (29) and (30), we see that they are almost 
adjoints of each another, with differences (relating to the ordering of the Tn and the position 
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of V) which are irrelevant to the statistics of C„ as n — >■ oo (since the %, are iid) . Noting that 
the decay rate A in (3) can be defined as the largest singular value of these operators, and 
that operators which are adjoints of each other have the same singular values, we conclude 
that the decay of the forward and inverse maps are identical. 

The prediction (9) for the decay rate in the locally controlled regime applies to flows that 
are sufficiently mixing. The precise notion of what sufficiently mixing means is not entirely 
clear, but we expect it to be close to the dynamical-system notion of exponential decay of 
correlations [35] . One property which plays an important role is that of transitivity of the 
stretching. This property, which is important for the large-deviation form (5)-(7) of the 
stretching statistics [36], requires that the tangent map An leave no deterministic directions 
invariant. As noted earlier, not all possible 3D extensions of the 2D sine map satisfy this 
property. To demonstrate this, we have considered the map 

Xn+l^ Xn+asm{yn + (l)n): Vn+l ^ Vn + b Sm{Xn+l + i>n) : = ^n + C sin(y„+i-F(^„) , (31) 

as an alternative to (14). The map (31) does not stretch exponentially in the vertical. 
Consequently the concentration field takes a highly anisotropic form, very different from 
that obtained for (14): the vertical gradients of concentration are much weaker than the 
horizontal ones, as illustrated by Fig. 10 which displays the concentration from a scalar- 
decay simulation using (31) and k = 10~^, and the scalar decay is controlled by horizontal 
stretching. In this case, the decay rate is obtained from (9) by considering the horizontal 
components of (31) only. 
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